FUNCTION medsigma2, inarr, sigma=sigma, dark=dark, gain=gain, average=average, im_stack=im_stack, out_dark=out_dark, help=help
;+
; NAME:
;       MEDSIGMA
; PURPOSE:
;       Compute the median at each pixel location across a set of 2-d images,
;       optionally ignoring outliers a specified # of sigmas away from the
;       median.
; EXPLANATION:
;       Each pixel in the output array contains the median of the 
;       corresponding pixels in the input arrays, with outliers removed
;       if desired.   Useful, for example to 
;       combine a stack of CCD images, while removing cosmic ray hits.
;       All input data types are converted to floating point.
; CALLING SEQUENCE:
;       outarr = MEDSIGMA2([inarr, sigma=sigma, dark=dark, gain=gain,
;                  im_stack=im_stack, out_dark=out_dark, /average, /help])
; INPUTS:
;       inarr  -- A three dimensional array containing the input arrays to 
;                 combine together.  Each of the input arrays must be two 
;                 dimensional and must have the same dimensions.  These arrays
;                 should then be stacked together into a single 3-D array,
;                 creating INARR.  INARR can also be a list of 2-D images
;                 from which the routine will create the the 3-D array.
;                 If INARR is not provided, the user is prompted to select
;                 a list of files to use.
;
; OPTIONAL INPUT:
;                
;       sigma  -- A sigma value; number of standard deviations away from median
;                 to accept.  Any one pixel in the set that has a value
;                 equal to or greater/smaller than the median+/-sigma,
;                 will be removed from the final median calculation. The
;                 meaning of sigma changes when a dark image is provided.
;
;       dark   -- If provided, this 2-D array will be subtracted from all the
;                 data files before the median is applied. The sigma then
;                 becomes the number of photon_noise/gain values above the 
;                 median to use to reject pixels
;                 (rejection_value = sigma * SQRT( (ADU-DARK)/gain )).
;
;		gain   -- The scaling factor between counts per pixel and the parent
;                 gaussian distribution; e.g., ADUs/electron for CCD images.
;                 Only used if DARK is set, to scale dark-subtracted counts
;                 before assuming std_dev=sqrt(counts).  Default = 1.0
;
;       average-- If this flag is set, the final image will be the AVERAGE of 
;                 all the included pixels instead of the MEDIAN.
;
;       help   -- If this flag is set, the calling sequence is displayed and
;                 the function returns immediately.
;
; OUTPUTS:
;       outarr -- A 2-D single-precision floating point array, with each 
;                 element equal to the median of all values at that location in
;                 the input image stack It will have dimensions equal to the
;                 first two dimensions of the input array.
;
; OPTIONAL OUPUT:
;
;       im_stack -- The original image stack can be returned via this
;                   keyword, e.g. if INARR is set equal to a list of file
;                   names.  (Otherwise, will be equal to INARR.)
;
;		out_dark -- The dark frame (if used) can be returned via this keyword,
;                   useful if DARK is set to a filename instead of an
;                   image.  (Otherwise, will be equal to DARK.)
;		
; RESTRICTIONS:
;		Requires IDL v5.6 or greater (uses DIMENSION keyword to MEDIAN).
;        Calculations are always in single-precision floats (since MEDIAN
;        always returns FLOAT data type)
;
; EXAMPLE:
;       Suppose one wants to combine three floating point 1024 x 1024 bias 
;       frames which have been read into the IDL variables im1,im2,im3
;
;       IDL> bigim = fltarr(1024,1024,3)        ;Create big array to hold images
;       IDL> bigim(0,0,0) = im1 & bigim(0,0,1) = im2 & bigim(0,0,2) = im2  
;       IDL> avgbias = medsigma(bigim, sigma=2)
;
;       The variable avgbias will be the desired 1024x 1024 float image
;       with outliers of at least 2 sigma ignored for median
; PROCEDURE:
;		Convert the input array to floating point.  If SIGMA is
;		given, set outlier values to NaN.
;		Apply MEDIAN over 3rd dimension, using the /EVEN keyword
;       (which averages the two middle values), since this avoids 
;       biasing the output for an even number of images.
;		Any values set to NAN (not a number) are ignored when computing the
;       median.    If all values for a pixel location are NaN, then the median
;       is also returned as NAN.
;
; MODIFICATION HISTORY:
;       Written by Michael R. Greason, STX, 12 June 1990.
;       Don't use MEDIAN function for even number of images.
;          W. Landsman Sep 1996
;       Mask added.  RS Hill, HSTX, 13 Mar. 1997
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Use /EVEN keyword to MEDIAN    W. Landsman  September 1997
;       Rearranged code for faster execution   W. Landsman January 1998
;       Faster execution for odd number of images   W. Landsman July 2000
;       V5.4 fix for change in SIZE() definition of undefined variable 
;                W. Landsman/E. Young   May 2001
;       Use MEDIAN(/DIMEN) for V5.6 or later   W. Landsman   November 2002
;       Use keyword_set() instead of ARG_present() to test for presence of mask
;           parameter  D. Hanish/W. Landsman   June 2003
;       Modified from Astron Library medarr function to add the sigma keyword
;           S. Beland December 2004
;		Vectorized for speed; all operations and outputs are floating-point,
;           since MEDIAN operates in FLOATs anyway.   Dropped MASK and
;           OUTPUT_MASK.  Introduced GAIN.  Corrected to ignore SIGMA
;           cutoff at locations where std_deviation=0
;           N. Cunningham   June 2005
;-

On_error,2
;Check parameters.
if keyword_set(help) then begin
	print, 'MEDSIGMA2 CALLING PROCEDURE:'
	print, 'outarr = medsigma2(inarr, sigma=sigma, dark=dark, gain=gain,'
	print, '        im_stack=im_stack, out_dark=out_dark, /average)'
	return,-1
endif

 ; check if inarr was provided
 if n_params() eq 0 then begin
	; get a list of files to process
	inarr = dialog_pickfile(/multi,filter='*.fits', $
		title='Series of images to median')
	if n_elements(inarr) lt 2 or inarr[0] eq '' then begin
		message, '2 or more 2-D images required in MEDSIGMA2'
		return,-1
	endif
 endif

 ; check if inarr is a string array
 if size(inarr,/TNAME) eq 'STRING' then begin
	inlist=inarr
	narr = n_elements(inarr)
	; get the size of the array from the first image
	img = mrdfits(inarr[0],0,h1,/unsigned)
	s = size(img)
	ncol = s[1]
	nrow = s[2]
	type = s[s[0]+1]
	tempimg = fltarr(ncol,nrow,narr,/NOZERO)
	tempimg[*,*,0] = img
	img=0b
	for fn=1,narr-1 do begin
	   tempimg[*,*,fn] = mrdfits(inarr[fn],0,h1,/unsigned)
	endfor
	inarr = temporary(tempimg)
 endif
 oparr=float(inarr)

 s = size(inarr)
 if s[0] NE 3 then $                    ; Input array size.
		message, 'Input array must have 3 dimensions'

;check/initialize sigma
 if n_elements(sigma) eq 0 then begin
	sigma=0
 endif else begin
	if sigma lt 0. then sigma=0
 endelse

;get array dimensions
 ncol = s[1]
 nrow = s[2]
 narr = s[3]
 type = s[s[0] + 1]
 help,inarr

 ;are we taking sigma cut?
 if sigma gt 0. then begin
	dodark=0b
	;check whether dark file was passed
	 if size(dark,/TNAME) ne 'UNDEFINED' then begin
		if size(dark,/TNAME) eq 'STRING' then begin
		   ; a filename was passed -> open it
		   dfile=dark
		   dark = mrdfits(dark,0,h,/unsigned)
		endif 
		sd=size(dark)
		if sd[1] ne ncol or sd[2] ne nrow then begin
		   print,'Dimensions of dark image is different from the data -> skip dark subtraction'
		   dodark=0b
		endif else begin
			oparr -= rebin(float(dark),ncol,nrow,narr)
			dodark=1b
			if n_elements(gain) eq 0 then gain=1.0 else gain=float(gain)
		endelse
	 endif
	medarr=median(oparr,dimension=3,/even)
	medarr3d=rebin(medarr,ncol,nrow,narr)
	;calculate SIMGA cut levels
	if dodark then begin
		;negative values in MEDARR will become NaN after SQRT
		nneg=long(total(medarr lt 0.))
		if nneg gt 0 then begin
			print, 'MEDSIGMA2:  Attempting SQRT of '+ $
				strtrim(nneg,2)+' negative values in MEDIAN(INARR-DARK)!'
			print,'Sigma cut will NOT be applied at these locations.'
		endif
		photon_noise=sqrt(medarr/gain)
		wzero=where(photon_noise eq 0.,wzcount)
		cutoff=rebin(photon_noise*sigma,ncol,nrow,narr)
	endif else begin
		tsqarr=total(oparr^2,3)
		totarr=total(oparr,3)
		sdev=sqrt(tsqarr/(float(narr-1))-1./float(narr*(narr-1))*totarr^2)
		cutoff=rebin(sdev*sigma,ncol,nrow,narr)
	endelse
	;check which pixels FAIL, not pass (NaN don't count as failing!)
	;don't apply sigma cut where cutoff value=0 (e.g. all values are same)
	pos = where(cutoff gt 0 and abs(oparr-medarr3d) ge cutoff, count)
	if count gt 0 then oparr[pos]=!values.f_nan
	rejected=long(count)
endif

;take median or average
if keyword_set(average) then outarr=avg(oparr,3,/nan) else $
	outarr=median(oparr,dimension=3,/even)

;return original image stack and dark frame if requested
 if arg_present(im_stack) then im_stack=inarr
 if arg_present(out_dark) and n_elements(dark) ne 0 then out_dark=dark
;Correct INARR and DARK so that variables aren't changed by this routine
 if size(inlist,/TNAME) ne 'UNDEFINED' then inarr=inlist	
 if size(dfile,/TNAME) ne 'UNDEFINED' then dark=dfile
 if sigma gt 0. then print,rejected,ncol*nrow*narr,format='("Rejected ",I,"  out of ",I," pixels")'

 return, outarr
 end
